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DESCRIPTION OF STATISTICAL METHODS AND A 
ROUTINE FOR DETERMINING THE PARAMETERS OF 
A ^30DEL IN PROCESSING EXPERIMENTAL RESULTS 

D. A. Usikov 

USSR Academy of Sciences, Institute of Space Research, Moscow 


The selection of the optimum parameters of a theoretical model and /2^* 
determination of the errors in them due to errors in physical measurements 
is an important stage in pjdbcessing an experiment. Besides, in processing 
experimental data it becomes necessary to evaluate the conformity of 
theory with the experiment. The routine described in this work solves 
these problems. The user wishing to process a specific experiment with 
its help has only t/ write a subroutine for calculating the function of 
the specific model. 

In compiling this routine attention was concentrated on assuring 
reliability, algorithmic speed and convenience. The routine extensively 
utilizes formatted printing and diagnosing possible errors in the input 
data. This paper describes in detail the specification sequence for 
the input data and the format of the calculation results. Necessary 
information on statistics is presented in a special chapter. / 

The programming language is FORTRAN, and the routine has been 

entered as a module in the routine library of the SOFI video display 

(■) 

processing complex. 

INTRODUCTION ^ 

i! 

II 

// 

The described routine is intended for processing ex|/eriments by 
the method of least squares or the maximum likelihood method. The 
class of selected functions is arbitrary, the required parameters may 
enter nonlinear ly. To use the routine it is necessary to program the 
function calculation block in each specific case. 

The result of the routine is a printout of a set of tables; 


♦Numbers in the margin indicate pagination in the foreign text 
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1. The optimum values of the parameters of the model; 

2. A covariance matrix of parameter errors; 

3. The theoretical curve optimally describing the experiment; 

4. The minimum statistical sum# i.e. , the sum of the squares of 
the differences betveen the experimental and theoretical values of the 

functions calculcafed for the best choice of parameters. 

1 1 ' ' 

j) 

Acctirding to the chi-square criterion, the value of the minimum 
statistical sum makes it possible to select competing models as well 
as to determine the degree of correspondence of the model and the 
experiment. 

Chapter I. ACCESS TO THE ROUTINE 
Sec. 1. Data Input 

The experiment is processed by the method of least squares. The 
sum of the errors (statistical sum) is minimized: 

- j I ^ 

h-t fe 

where M is the number of experiments; is the coordinate of the 
i-th measurement, for example, the instant, length, etc.; f_(Z.) is 
the experimentally obtained values at points Z^; f^(Z^j^,x) is the 
theoretically predicted value at point Z^; X is the vector of the /4 

selected parameters of theory. We denote the dimension of vector X 
as N(N = dimX) . is the error of the i-th experiment (one standard 
error) . 

The routine looks for the values of X at which the statistical sum 
(1) is minimal. The accuracy with which the minimum is sought is given 
by a special parameter related to the statistical nature of the problem. 
Determination of the parameter is described in Chapterll. It is 
al^sumed that the statistical sum (1) has one minimum. If there are 
several minima, a local minimum is determined, depending on the initial 
approximation of the parameters. 
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In the description of the quantities that follows, the notation 
employed corresponds to the identifiers in the routine. The quantities 
are entered in the same sequence as they are described below. 

Quantities Entered Into the Main Routine ; 

N - dimension of the space of the X parameters. Format 16. Restric- 
tion, N 42 O. 

IPE - printout label. IPE = 0 - normal printout mode. IPE = 1 - test 
printout Tode. Format 16. 

Quantities Entered by ENEXPE Subroutine ; 

M - number of test points. Format 16. Restriction, M<400. 

Z (M) - array of M numbers - measurement coordinates . Format 5ElS . 7 . 
EXPE(M) - in notation of (il) - fg(Z^) - array of experimental values. 
Format 5E16-7. 

VEXPE(M) - in notation of (1) ^ ~ array of experimental errors. 

Format 5 E16.7. 

Quantities Entered by SPQINT Subroutine ; 

X(N) - array of parameter input values,, The closer the input parameters 
approach the optimum values the faster, in general, is the minimum 
found. Format 5E16.7. 

AVD(N) - array of input values of parameter errors. These quantities 
are required by the routine for the initial selection of the gate 
circuit on which the first and second derivatives are calculated. 

The quantities are subsequently modified by the Routine as the Work 
proceeds. It is recommended to take AVD(i) ^ O.lX(i), i.e., take 
the errors at approximately 10% of the input values of the respective 
parameters. 

Attention: AVD(l) should not be taken equal to zero! Format 5E16.7. 
MARK(N) - qualifying array. If MARKi(i) = 0 the given parameter varies. 

If MARK(i) = 1 the given parameter is reinforced and taken equal to 
the input value of X(i) . Format 7211. 


Quanta, ties Entered Into the Main Routine: 


EPS - accuracy of determination of the minimum, EPS is usually taken 
equal to 0.1. The precise definition of EPS is given in Chapter II. 
Format E16.7. 'I 

Sec. 2,, Printout of Results 


Below is presented ti^e printout sequence of results for the case 
IPE = 0, i.e., when a test is not envisaged (sefe Sec. 1). 

We shall illustrate the routine printoit with a concrete example. 

The heading is printed indicating the model employed. This is 
followed by the input data: the EXPE(M) experiment array, the array 
of test errors VEXPE (M) , the coordinates of the test points Z (M) : 
EXPERIMENT APPROXIMATED BY I ORDER POLYNOMIAL 


iPf-V - TEST PRINTOUT YES 


POINT NUMBER 
1 
2 
3 

M 

5 

6 
? 

8 

9 . 

10 


experimental value 

0l 

01 

0.3A09999S 01 
0.Sf9r00O* 01 
0*5ft0999f| 01 
C.SfVftOOOf 01 

0. ?99/^oOOf 01 
G.9nfl999V| 01 

1 . - 


experiment error 

0.5999958f»*02 
0 • 999999e8'»02 
0.9999998|9«2 
0.9999998EB02 
0 •9999998E«02 
0.99999988^02 
0.9999«98E!so2 

0 ‘9999998EBQ2 

0 . 9999998EB02 

0,999999881102 


COORDINATE 

o.iooooeof 01 
o.aoooooof el 
o.3ooooeo| 8l 
o,«oooooo| 01 
0 'Sooooool 01 
0,6000000| 0l 
o^?v^oeooo| 01 
0 , 8 o 00 o 00 i 0 l 

0 . 9000000| ol 
O.ldOOoOOi 08 


This ig followed by the input values of the parameters X(i) and 


the approximate errors of the parameters AVD(i) 


INPUT PARAMETER 
0 . 0 
0.0 


APPROXIMATE PARAMETER ERROR 

e.iooooooE 01 

0.1000000E 01 
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Some parameters may not v^'.ry (which is indicated by units in 
the MARK(i) array), and the routine therefore transfers to the internal 
list, which gives only the varied parameters. Correspondence between 
external and interna numeration is indicated in the following table: 


DATA PRINTOUT EVERYWHERE ACCORDING TO INTERNAL NUMERATION 
CORRESPONDENCE OF PARAMETERS ‘ 


INTERNAL PARAMETERS 
-2 


EXTERNAL PARAMETERS 
r 


Next the EPS value is stated: 

accuracy of DETERMINATION OF MINIMUM 
-scale ONE CHI-SQUARE DISPERSION 


With this the input data printout ends and the routine transfers 
to calculation; 


1 

2 

3 


THEORETICAL CURVE DEMONSTRATION 
EXPERIMENTAL VALUE ” 

■ ... T.S, ..... 






POINT Y 

^ I 


f* • 


• «» W •»«»«« 


i 

X 

I 

X 

X 

I. 

I 

t 

r 

t 


FIRST INPUT 

PARM'TrR=«‘ 

1 

2 





I 

t 

I 

I 

I 

J 

? 

I 

I 

J 


0 . 1 ft 1 0 i 
0 . i f 90S 
o.3(noe 
0 . 3^901 
0 .5010 1 
0 . 5990 | 
0 . 701 OE 
0.7990E 
O.9010| 
O.9990B 


01 

Cl 

01 

01 

01 

01 

01 

01 

01 

01 


’ 0 * 10 ". o| 

* 0.19906 
' 0 . 3010 E 
0.39901 
0 .soioj 

0 * 5990 f 
0.701 Of 
0 . 7990 E 
0 ‘ 9010 E 
0 . 9990 E 


0 «» • 
^ ft 

I 

03 i 
03 X 


03 

03 

03 

03 


03 I 


•r mtm m ^ f* m mm irnpUr fVmmnrnmmff 


7 • ‘STATISTICAL SUM 

ETER VALUE LAST BASE .... 

.0 0,30120076 00 


ORi’GflVAL PAGE |: 
OF POOR QUALITY 



0,4aSft292i-01 


As the statistical sum decreases the pro^jrram prints out 
intermediate results ; 


•INPUT FOLLaWINa'NEWTON METHOD STEP 
o.*«*9‘«7e 01 - STATISTICAL SUM 
PARAH^T€R PVAWW.likilE 

2 0.99934226 00 


0.77647262-04 


ORIGINAL rAC- 
OF POOR 


1? ^ NUMBER OF ACCESSES TOSTATSUM BLOCK 
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When the optimum value the parameters has been determined 
the table of theoretical and experimental values is printed out| after 
which the routine gives the minimum statistical sum, the Fisher 
matrices, and the covariance and correlation matrices; 


t POINT' t 


I 

I 

I 

I 

l 

I 

I 

I 

i 

I 

<•»«»»« « 


THEORETICAL CURVE DEMONSTRATION 
EXPERIMENTAL VALUE 
COMPARISON WITH EXPERIMENT 
(THEORY MINUS EXPERIMENT DIVIDED 

1 J 

n 


1 

i 

S 




Ull!E.DlUD|ejY-li§lAWOR 



0. loose 
o.zftozi 
o*s<'OEi 

o*4ftoie 

0- Sf'Oftf 
0 » 6 noft| 

0*6999f 

e,799Re 

Q.8998e 

o.999?e 


01 

01 

01 

01 

01 

01 

01 

01 

01 

01 


I 

I 

I 

t 

1 

! 

I 

! 

I 

I 




0,10101 

0.1 990 I 

0‘3010| 
0,39906 
0*50106 
0,59901 
0» 70106 
0,79901 
0.90101 
0. 99906 


01 

01 

01 

01 

01 

01 

01 

01 

01 

01 


I 

1 

I 

I 

I 

I 

I 

I 


0.72736 

mill 

0 .96961 

riini 

0*1212E 

0,72736 


00 

01 

00 

u 

00 

01 

00 

01 

00 


I 

I 

I 

I 

t 

I 

I 

I 

I 

! 




MINIMUM SEARCH END 

0 .i, 8 ^. 82 UE M - STATISTICAL SUM 
PARAMETER PARAMJJIR^VALUE^ - 

2 0.f9f342?E 60 


4\SJ7M^6-6S 

ft. 776^7266*04 




NUMBER OF ACCESSES TO STATSUM BLOCK 


FISHER MATRIX 

1 1« e,927t255E u5 

41?“ 0.51002791 «6 

“ 0,35586966 c7 


8 

8 

8 


. COVARIANCE MATRIX • 

1 I" 0 . 5079759E»,.ft4 


CORRELATION MATRIX 
K 1 1“ e.iortO''oo| Cl 

K 1 ?» "0,88759626 00 | 

K 2 7“ 0,1060000E ,1 " . 

The routine cycle is complete. When the calculrihion ends control 
is transferred to the start, and the routine requires input of a new 
set of data. 

Sec. 3. Function Calculation Subroutine 


The function calculat^^dh subroutine^ FUN(i, Zi, N, X, FT) is the 



only subroutine which has to be changed as applied to the solution 
of a specific problem. The subroutine provides for calculating the 
function value for the given parameter values of model X(N) at point 
Z^. For example » for a linear model the function is calculated a/;jcor- 

ding to the formula FT = X(l) + X(2)xZi + X(3)xZi^2 + The PUN 

subroutine operates in two modes. The respective mode specifies the 
value of the identifier i. At i - 0, the subroutine does not perform 
calculations and only prints out the heading 

EXPERIMENT APPROXIMATED BY 1 ORDER POLYNOMIAL 

At i = 1, calculation of the function takes place. 

The formal parameter N in the subroutine access specifies the 
number of parameters of the model (varying or nonvarying in sum) * 

The parameters i, Zi, N, Z enter the FUN subroutine from the main 
routine, and the value of the function FT is transmitted back. 

At points lying far from the minimum the search mode may generate 
values of the X parameters that are unusable in FUN subkoutine calcu- 
lations. To avoid such a situation a special message contingency is 
provided for. Before starting FUN subroutine calculations it is 
necessary to check xvhether the incoming values of the X parameter 
are permissible. If they have gone beyond the limit of permissible 

18 

values the calculation is cancelled and FT is assigned the value 10 

18 

When the pilot routine receives PT equal to 10 it takes steps to 
enter the domain of permissible values of the X parameters. 

There are different ways of determining the permissibility of 
the X parameters. 1) X is verified prior to the calculation by means 
of a test consisting of a set of inequalities. 2) The permissibility 
of X is verified during the calculations. 3) A special FORTRAN device 
is employed: the possibility of transferring control to a specified 
place in the subroutine. 






Sec. 4. Procedure For RecalXinq Subroutine from SOFI Library 


The subroutines are cataloged in the object module library of 
the SOFI complex (1), The routine is generated in the following way: 

/f JOB NAME 

// PAUSE ASSGN SYSRLB,X' 190 ’ disk 102 
// OPTION LINK 
// EXEC FFORTRAN 

CALL MODNEU 
STOP 
, END 

SUBROUTINE FUN(I, ZI,N,X,FT) /lO 

DIMENSION X(20) 

CAI/L FUN. . . (I,ZI,N,X,FT) 

RETURN 

END 

SUBROUTINE .^tiN . . . ( I , Z I , N , X , FT ) 

— 

subroutine bo<|y 

/* 

// exec lnkedt 
// exec 


input data 


The MODNEU subroutine contains the body of the main routine (texts 
of the routine in FORTRAN are given in the Appendix) . In FUN are 
cataloged subroutines far calculating different functions. At present, 
the subroutines FUNPOL and FUNEXP, which calculate polynomials of 
the Nth order and the sum of exponents, respectively, have been cataloged. 
A special FUN subroutine is written for each of the calculated sub- 
routines. The FUN subroutine transfers control from general access 
of the main MODNEU routine to the FUN function to the concrete FUN... 

# 

subroutine. 




If a calculation of the same type is performed for many variants 
according to some PUN . . . subroutine , an absolute model can be generated 
and put into the SOFI complex according to the procedures cited in 
work cu- Thus, for example, at number 21000 in the SOPI complex there 
is located a routine with a PUNPOL subroutine, and at number 21001 in 
the absolute library of the SOPI complex is a routine with a FUNEXP 
subroutine . 


Chapter II . ACCURACY OF DETERMINATION OP MODEL PARAMETERS 


Sec. 5. Normal Limit of Likelihood Functions 

IIP [■ ■ C I I n il 1 11 I II i n[ i h im iii m i w i . h i. ;i i i i jiu.i ip h 


In this chapter are described statistical methods of processing 
experiments involving normal approximations of likelihood functions. 

The domain of applicability of normal descriptions of an experiment 
is considered at the end of the chapter. 

The likelihood function 1 (T/E) of a model T with respect to a 
given experiment E is determined according to Bayes equation; 

'i\ 

/P(Bn)P(T)cfT ( 2 ) 


Here, P(E/T) is the probability density of the realization of the 
experiment E, provided the parameters of the model are T; P(T) is 
the a priori probability density of the parameters of the model. 

The likelihood function is sometimes called the "a posteriori" proba- 
bility density of the model parameters, and P(T) is the "a priori" 
density. 

It is assumed that for sufficiently representative experiments 
variations of the a priori density P(T) are insignificant in comparison 
with the peak of the likelihood function at the maximum point 1(T/E). 

In these assumptions it is proved that, as the number of experiments 
increases, the function 1(T/E) tends towards a normal distribution 
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P, p. 2173 . A multidimensional normal distribution has the form: 


PfXj 




'f/i 


e^p(-j (x-Xcf ^ 


(3) 


where IM is the determinant of A; X*AX » b •* dimX is the 

dimension of the parameter space. 


The loatrix A is called the '’information” matrix (or the Fisher 

*-l 

matrix) of the experiment; its inverse B = covX * A is the covariance 
matrix. Matrices A and B are symmetrical and positively determinate, 
i.e., X^AX>0 for all X, 




The normal distribution is given by two parameters; the mean value 


Xq ~/xP(X)dX 


(4) 


and the covariance matrix 


B = covX = - Xq) (X - XQ)^P(X)dX. 


(5) 


It is thus assumed that the experiment is fully defined if Xq 
and B are known. Calculation of Xq and B according to equations (4) 
and (5) is inconvenient algorithmically, and differential methods are 
usually employed. The mean values of Xq are usually found from the 
maximum condition of the likelihood function: 


PX; 


O. 


( 6 ) 


The Fisher matrix is obtained as the second derivative of the logarithm 
of the likelihood function; 


/i; 






ii) 


It is not hard to see that definitions (4), (5) and (6), (7) are 
equivalent. Determination of the parameters Xq and B (or A) from 
equations (6) and (7) naturally^ Jn volves numerical methods of looking 


10 


for the maximum of the likelihood function. This relationship will be 
investigated in greater detail in Chapter III in describing the 
algorithm. 

Description of an experiment by stating the parameters Xq and B 
is called "normal" description of the experiment. Normal description 
forms a sufficient statistic for any linear combination of input 
parameters. Let us determine 


KX, 


( 8 ) 


where X and Y are vectors, dimY = r, dimX = n, and K is the matrix 
of the dimension coefficient rxn. If vector X is distributed 
normally with the parameters Xq, B, then Y is also distributed normal- 
ly with the parameters: 


/13 


Yq = KXq, covY = KBK^. 


(9) 


The relationships (9) are also frequently employed in the case 
of a nonlinear relationship Y = f(X), the matrix K being determined 
as the factors of the linear term in a Taylor expansion of the function 
f(X) at point Xq: 

9j// j 

Sec. 6. Distribution of the Statistical Sum. The Chi-Square Criterion 


The most common case is when the errors of an experiment are 
distributed according to a normal law. The function P(E/T) (see 
equation (2)) has the form: 

P(E/T) = c exp(-3sfg - f^(X))^ £(fg - f^(X)). 


( 11 ) 


Here, c is the normalization constant; $3 is the covariance matrix 
of the errors of the experiment; 

' is a vector compounded of experimentally measured 
values; 


£ ^ 



OF POOR 
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/xj’k is a vector compounded of the theoretical values 
^ of the function in the i-th experiment, when the 
values of the theoretical parameters are X; 

M is the number of experimental points. 


(V '"/i 


Note that the X parameters in the model f^(X) may enter nonlinearly.. 
The function 




('f, 7/y 


( 12 ) 1}A 


is called the ’’statistical sum" of the experiment. Function (12) 
does not differ essentially from the function (1) introduced before. 
Function (12) generalizes (1) for the case of dependent experiments. 

Expanding S(X) in a Taylor series of X in the neighborhood of Xq 
and neglecting terms higher than the second order of smallness, we 
arrive at the normal description of the experiment; 

S(X) = S(Xq) + (X - Xq)'^P'S(Xq) + Js(X - Xq)^A(Xq) (X - X^) , (13) 


Here, ^S(Xq) is the gradient vector; 





A is the matrix of second derivatives; Aj^. = ^x.^x. 

It is assumed that S(X) has one extremum. Point Xq is found 
from the condition 


.?S(Xq) = 0. 


(14) 


In particular, if the parameters enter the model linearly, then 


f,j,(X) = FX, 

where F is the dimension matrix dimf x dimX (M x N) . 
the expansion (13) is exact, and 


(15) 

In a linear model 
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A = 


(16) 


The linearization procedure makes it possible to formulate 
a very convenient accuracy criterion for numerical methods of solving 
sets of equations. Suppose we have to solve a set of equations 


» X 2 f * • • » Xjj) — 0; i — 1^2,. ..^N^ 


(17) 


/15 


or in vector notation, F(X)X = 0. 


Specific requirements are imposed on the accuracy of determination 
of the parameters, namely. 


2 


/ X/ -X/af 

S'y 


< £ 


(18) 


i.e., the approximate solutions of X should differ from the exact X 
by no more than a specified quantity characterized by the given error. 
More precisely, condition (18) is equivalent to the requirement that 
the deviations of the approximate values must lie within the ellipsoid 
given by equation (18) . The parameter » states the degree of approxi- 
mation to the exact solution. If it is necessary to take into account 
the paired relationships of the accuracies, then condition (18) is 
replace by the condition 


(X - Xq)^X~'(X - Xq) , 

where is is a positively given definite matrix, 


(19) 


X = 


(O' 


^iO = 


\Xjjo/ 


The solution of equation (17) is apparently equivalent to the 
solution of the problem of finding the minimum, 


^ 2 T 

min >4Fj|^(X) " min F F, 


(21) 


This procedure, which is extensively employed in numerical methods, 
makes it possible to solve the initial problem by developed methods of 

13 


looking for the extremum of a function of many variables. However, 
if i; 

( 22 ) 

is adopted ,as a cohdition for attaining the required accuracy, the 
values of X which satisfy (22) will not, generally speaking, satisfy 
the inequality (19) . Let us now show that, proceeding from the 
inequality (19), we can formulate an equivalent inequality in terms 
of F. To the accuracy of the lineai; term, we have 

F(X) = P(Xq) + K(X - X),, (23) 


where 


K 


- 
ij 


X=X, 


By definition, F(Xq) = o» hence (23) involves only the 


linear 


term 


F(X) - K(X - Xq), 


(24) 


If solution (17) is unique, then K can be inverted, and there 
exists a matrix K Substituting X - Xq ~ K ^F into (19) , we obtain 




or 


F*^ (K s K^) .%• I (25) 

Matrix is positively definite, because £ is positively 

definite. Matrix (KIc k'^)''^ is also positively definite. Finally, 
the solution of equation (17) is equivalent i:o the solution of the 
problem in finding the minimum: 

min F^(Ki: k'^)“V, (26) 

X , 

and the condition of attaining the minimum (X - X.)^2~\x - X.) is 

T T —1 ^ ^ ^ 

equivalent to the condition F (KS K ) ps» 


Let us return to an examination of the properties of the 
statistical sum (12) . Since it is impossible to find such an Xq 
that (14) would vanish by numerical methods, it is necessary to 
formulate a criterion that would characterize the degree of approx- 
imation to the exact solution of the extremum problem. It is natural 
to require that the numerical method should be the more precise the 
higher the accuracy of the experiment. In other words, the accuracy 
of the search for the extremum should be related to the covariance 
error matrix. We recall that is the evaluation of the parameters 
of the model. The covariance matrix B = expresses the probability 
of deviation of the true vaUie of the parameters from the value 

of Xq obtained from the condition of the likelihood function maximum. 
The main statistical criteria are linked with the distribution 
instants of the statistical sum (12) . It is not hard to show that 
the statistical sum possesses a chi-square distribution, provided 


f(X) is linearly dependent on X. The mean value S =; 
central instant? 


M-N 


the second 








o 


At M - N of the order of lO or more, the chi-square distributiou 
close to its maximum point can be considered close to normal. Therefore 

if |s(Xq) - ^ ^ then the chi-square criterion is satisfied 

to a 68% probability. With the help of this criterion we can judge 
of the correspondence of theory and experiment, for example, corres- 
pondence of the experiment errors aT^ipdicated by the experimenter 
to the true experiment errors. For examp,Le, if it seems certain that 
the experiment is correctly described theoretically, while as a result 
of the search for Xg it has been found that S(Xg)< S -a , this can be 
interpreted as indicating that the experiment errors have been assumed 
too high. 

According to the chi-square criterion, the solution of equation 
(14) is considered satisfactory if 


Js(Xq - X)^a"^(Xq - X)« *if. 


(27) 


or. 


U 
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> 

I '■ 

/V . , / 

Here, X is the approximate value/ of Xq, * is a certain constant 
expressing the degree of confidence* At # s X, as noted before, /18 

the confidence is 68%. At i> 1, the confidence is accordingly 
higher. Usually • is accepted equal to 0,1, 

Since Xq is not known, criterion (27) is not constructive but, 
taking into account that S = A(X - Xg) , we can obtain from (27) 
an easily computable criterion: 

% r S*^a"^ VS . (28) 

I 

The gradient V S is computed at point X in the search for X^ analy- 
tically, if the explicit^ form f,j,(X) is known, but more often approxi- 
mately. The matrix of second derivatives is also calculated in the 
course of the computations. Criterion (28) is especially convenient 
because calculation of r S and A are essential in the Newton method 
described further on. ' 

Sec. 7. Condition for Ending the Minimization Proces s 

The experimenter often does not know the absolute values of 
the experimental errors iT^. All he knows is their relative course 
from experiment to experiment. Consequently, the quantity <r in (28) 
is indeterminate. In that case we can make use of the fact that 
* = and, dividing both parts of (28) byV S , make use of the 

criterion . 

h ^ s^^a"^ vs^ Vs. (29) 

The unknown quantity S in (29) is substituted by the current 
statistical sum S(X): 

JstS^A-lf S« ^ ijs. (30) 

Criterion (30) is, obviously, weaker than (29) at S(X) 5J>S, However, 
in practice calculations show that in the overwhelming number of cases 
the quantity Jj V S A ■^S when it is far away from Xq increases faster 
than ^S, and criterion (30) does not lead to false stopping far 
from the extremum point. The criterion of stopping the minimization 
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process (30) is employed in the program described further on. 

I; 

Sec. 8. Statistical Description of Parameters of "Poorly Con^ /I 9 

ditioned" Models and Experiment s 

There are at least two cases when normal description of an 
experiment is unsatisfactory. The first, which is frequently en- 
countered in practice, is when it is impossible to reconstruct the 
parameters of the model from the experiment owing to low accuracy 
or inadequate statistic of the experiment. Such a situation arises, 
for example, in attempting to determine a great number of parameters 
of the model from a small number of experiments. The problem of 
describing the experiment that arises in this case will be discussed 
in more detail further on, , 

The second case, which is in a sense diametrically opposite to 
the first, is encounteredi. in processing a great number of highly 
accurate experiments. In this case it is usually found that either 
the model does not adequately describe the experiment or that the 
experiment errors are stated imprecisely. For example, the correlation 
between individual measureigents acquire major significance, or small 
systematic experimental errors become decisive, i.e. , the "trifles" 

/'■' '■ v' 

usually ignored/ but which in rich statistics restrict the attainable 
accuracy of reconstruction of parameters. This case can be detected 
from the chi-square criterion. When the models of theory and expe- 
riment differ from the real-life experiment the minimum of the 
statistical sum differs substantially from the theoretical value, 
which can be a "trouble" indicator. The situation is described in 
detail in work [3j. u 

We shall assume that the model is accurate and the errors of 
the experiment are given correctly. Let us consider the effect of 
the nonlinearity of function f(X), According to Bayes' approach [2], 
confidence, that the .theoretical parameters lie within the X domain 
is found from the formula 

D(X) = j/l(T/E)dT, 

X 


(31) /20 




« 


where 1(T/E) is the likelihood function of the experiment (see 
equation (2)). 


If the likelihood function is normal, then to calculate the 
integral (31) over any domain X, it is sufficient to know the 
normal distribution parameters: the vector b means and the covariance 
matrix. These parameters form a sufficient statistic. But if the 
likelihood function is not expressed by a normal distribution a 
situation arises which is conveniently illustrated with the help f 
of curves of the levels of the likelihood function: 



and X 2 are parameters of the model. Suppose that the likelihood 
function has only one maximum, and ^10 and X 2 Q are the parameters at 
which it attains that maximum. The curves show the solutions 
KX 1 X 2 ) = const for various constants. The level lines always form 
ellipses in the neighborhood of the maximum, but farther away from 
the maximum point, when the model f(X) is hoalinear, they are no longer 
ellipses. It is not hard to show that, for each confidence D there 
is such a d that in integrating over the X domain defined by the 
condition 

l(X)Sd, (32) 


we obtain 



-0 = ^(X)dX. 

Usually some meaningful confidence is assigned, for example, 0.68 /21 

or 0.99 


Let us denote a normal likelihood function with parameters Xq 
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B as 1 (X ) t The fundamental thesis of normal description of 
%ii;)experiment is formulated as follows. If 



-/= 

•'x 


?(X)dX 




(33) 


where D is a significant confidence, and the integration domain X 
is found according to (32), normal description of the experiment is 
assumed satisfactory. 


In practice condition (33) can be obtained by the Monte Carlo 
method. The integral ^l(X)dX in (33) can be determined according to 


the formula 


^l(X)dX = l(X)dX. 


(34) 


X^ is played according to the density of 1(X) (an algorithm for 
modelling abnormal distribution is described in wdrk [d] ) , The 
quantities , which are an evaluation of the integral (34) , 

JL Aa 

are are introauced into the summator. 

Chapter III. DESCRIPTION OF THE ALGORITHM 

Sec. 9. The Modified Newton Method 

The main task of the experiment processing algorithm is to find 
the minimum of the statistical sum S(X) (12). The modified Newton 
method [4, 5, 6] used for this is based on local quadratic interpola- 
tion. Suppose the statistical sum depends upon the parameters in 
the following way; 

V; 

S(X) = Sq + - Xq) + %(x - Xq)^A(X -■ . (35) 

Here, Xq is the only minimum point, and ^the gradient VS is taken at 
point X«, i.e., it is zero; H 


S(X) - Sq + ?5(X - Xq)%1x - Xq) . 


(36) 


/22 
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We take the gradient of both sides of (36) ; 

V S(X) « A(X - Xq) . (37) 

We multiply both sides of (37) by A and obtain the formula of 
the Newton algorithm: 

Xq=X-a”^>S. (38) 

If the representation S(X) (35) is exact, then the minimum point Xq 
is found, starting from point X, in one step according to equation 
(38) . Besides the Newton Method, the minimum can be found in a finite 
number of steps by means of the algorithm of the conjugated gradient 
method [7, 8] . The simple method of gradients, or quickest descent, 
does not generally speaking, converge at Xq in a finite number of 
steps [4]. 

If the expansion (35) is not precise, we may find that the function 
at the new point after a step by the Newton method is greater than 
in the preceding point: 

S(X - a"^ vs) > S(X). (39) 

To avoid such a situation and assure that the minimization process 
yields a monotonous decrease of function S(X), equation (38) is 
modified: 


Xq = X - «. a"^ vs. (40) 

If the simple Newton method fails to work on some step, « is selected 
in such a way that 

S(X - 4 a"^V’S) SS(X), (41) 

This device is known as the "modified" Newton method. 

Sec. 10. Internal Scaling / 

It is well known that a programnier, besides Selecting a good 
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algorithm, shoiild take account of the specific features of computers. 
Algorithms are reliable and sufficiently universal only when account 
is taken of cases going beyond the machine's digital grid. One of the 
methods of avoiding overflow during calculation is described in Sec. 3. 
Here we shall exantne scaling as related to the fact 'that the theory 
parameters X can be expressed in arbitrary physical units. 


The scale adopted for each variable X^^ is the evaluation of the 
error, or fiore precisely, the quantity; 




V 


dx, 


(42) 


This device makes it possible to make all the variables X^ dimension- 
less. The more precisely the quantity X^ is known the larger its 
value in the dimensionless form (42) . The scales of ^i are verified 
in each step of the Newton method. The transfer to the true measure- 
ment units is carried but at the end of the calculation, that i^, after 
the minimum has been found. 


Sec. 11. Block Diagram of the Routine 
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The routine and subroutine texts are presented in the Appendix. 
The synunetric matrix inversion subroutine was taken from the col- 
lection of algorithms [9]» algorithm N. 66b, and translated into 
FORTRAN. 

Sec. 12. Algorithm of Approximate Calculation of Gradients and matriit 
of Second Derivatives 


To reduce the user's preparatory work to the minimum, the routine 
is so devised that it does not require special programming of the 
first and second derivatives of the function with respect to the 
parameters of the model. These quantities are found in the FISHER 
subroutine by forming finite differences, f 1 /25 

t 

The following formula is used to compute the gradient: 
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if 

0x,- 


s{)(; 1‘A x; ; - s(xi - 


2 AX,- 

The second derivatives are calculated this way; 


'-<^J y 9x;7Xj 


9^s _ sf%: ^AX;,xj * .q “j) 


AX; aX, 






7‘s 


S( XI -aX;J -2<S(X;) f S(X; ~ AX/ J 

(AX.-f 


The problem of selecting the base A X. is resolved by taking 
as ^ X^ the evaluations of the parameter errors (42). 

Sec. 13. The Algorithm of One-Dimensional Minimization 

In passing from calculations by the Newton method to the modified 
method it is necessary to determine the minimum of function S(X - « A ^vs) 
as a function of the parameter • . The algorithm of the corresponding 
one-dimensional minimization is constructed as follows. By the time 
of reference to subroutine ODM the values of the function S(«) are known 
at two points a.; 


1) S2( «2 . 0) = S(X ) ; 


>-l. 


2) SgC .3 = 1) = S(X - A -^VS) . 


If = -1) < S( ) • a new point 

• = - (a„ - •, ) 2 


(43) 


is chosen. If S ( a ) < («j^ ), the points are redesignated, which 

cahn be written using ALGOL notation : 


“2 == • 1 ' *1 ' 


and theroutine transfers to executing the operator (43). 


This occurs until S(« )>S(*x)* redosignation occurs: 

di I J* <df \ df d \ 

Si ; A S, 

The aggregate of points «2,*3 such tat S(«^) >S(«2'» S(«2)> S(«2) 
is called '’canonical”. 

After obtaining the canonical triplet cf numbers the algorithm 
transfers to looking for the extremum by the quadratic interpolation 
method according to three points. The minimum point a is determined 
from the formula 


iff if ^3 (t 

d — * "'* 

et, - Ha 

whei:e 

Ct, = Sr (dj '•djJ ; 

O2 * Si (d, " dij y 
Oj * Si (d. 

t, -(da *d3 )/^ ' 

* ( 1 ^, ^dj)/ ^ ; 

f,^(d,td 3 )/ 2 . 


It is not hard to show that point cv. always lies within the 
interval(*2^ < From the aggegate of four points « , "l ' “ 2 3 ' 

S, S2, Sg, three such points#^ ®2' selected 

so that the requirement •]<«2< *3, S^Sj, The cancxiical triplet 

thus obtained is used to conpute the next aK>roximation of a according to 
the procedure described above. The condition for halting ibe one-dimensional 
minimization process is a one-dimensional case of the general criterion described 
in Sec. 7. 


1 


The described algorithm is highly effective. As numerous /27 

I 

calculations show, the average number of selections of canonical I 

, , 

triplets of points is about five. | 

Sec. 14. Drawback.i|\ of the Algorithm of the Modified Newton Method. f 

The Gradient Line Method. ■ 

5 

i 

Formula x,. = X ■* a“^V S, on which the Newton method is based, | 

is obtained from the condition VS(X«)~ 0. However, the condition is I 

. ^ I 

satisfied not only by the minimum points, but by the saddle points as i 

well. If in the iteration process point X approaches such a saddle | 

point the direction a”^Vs will x>oint to the saddle point, IK some | 

configurations of point X, going away from the saddle point with the | 

help of the one-dimensional minimization process may require a lot j 

of calculations. Let us explain this with the example of two variables, | 

Xj^ and X 2 . Suppose X = 0 is a saddle point: | 

I 

I 

i 


! 


The lines of level S (X) = const are marked with arrows indicating 
the direction of increase of function S (X) . If point X is located in 
the octants (X 2 > 0; Xj^ < 0) or {X 2 C 0 ; Xj^> 0) , the direction of 
one-dimensional minimization is towards point X = 0 , and the X itera- 
tion will remain close to the saddle point as long as the next point 
does not move to octants 0 ; X2*<0) or 0; X2'^®^ es a 

result of errors in the computation. The greater the number of 
variables the higher the probability of X approaching the domain 
of the saddle point. To get out of the saddle point the routine 
must carry out many iterations , which is precisely the principal 
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drawback of the modified Netwon method. Work |j5][ suggests ways of 
preventing X from occurring in the saddle point domain. The idea 
of the method is for matrix A, which is not positively determinate 
in the saddle point, to be replaced by a positively determinate 
matrix, for example. A* == (A^A)^. Matrix A can also be reduced to 
diagonal form by similarity transformation, and all negative elements 
be substituted by positive ones. Work [6^1 presents the results of 
test computations which show that substantial acceleration of the 
operation of the algorithm can be achieved in this way. It is worth 
noting in passing that an algorithm in which the statistical sum is 
in the linear approximation (16) always has a positively definite 
matrix. This, apparently, explains the relatively high effectiveness 
of this method. 

The second shortcoming ot} the modified Newton method is due to 
the need to invert the A matrix. If two parameters of a theory cor- 
relate strongly, or if the number of experiments is smaller than 
the number of parameters of the theory, inversion of the A matrix may 
prove difficult. In general, the following device can be employed; 
instead of the A matrix obtain an allied matrix of A*, i.e,, a matrix 
made of algebraic complements A^j^, and compute the determinant A 
separately. The structure of the allied matrix is; 


hi 

hi 


hz 

CVJ 

CM 

... A ^2 

hn 

c 

CM 

... 


An allied matrix always exists. It is related to the inverse matrix 
by the formula A* = A~^^A^. The direction of A~^Vs will, apparently 
coincide with the direction of A*Vs. 
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The two mentioned shortcomings can be overcome by further modifying 
the Newton method. In particular, in the following version of the 
described routine it is suggested to nse a new algorithm of minimi- 
zation along the gradient line, which is also based on local quadratic 
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approximation of the function. The idea of the algorithm is to 
proceed from the given point X along the gradient line and thus 
find the one-dimensional minimum. Moving along the gradient line, 
unlike moving in the direction of the gradient, makes it possible 
to attain the minimum in a finite number of arithmetical operations. 
The path of motion is described by the equation: 


^ X _ 

at ~ 


(44) 


for the initial condition X(0) “ ^in* ^ parameter. Since 

'iT S = A(X - Xq) , from (44) we obtain the equation 

p=A(X-X„), S i'. (45) 

which has the solution: 

X(t) = + (e^^ - I)a“^VS(X^j^) . (46) 


At t = 0, X(0) = ^ matrix is positively definite, 

then X(-oo) = - A ^V's, which coincides with expression (38) 

for Xq obtained by the Newton method. Thus, the parameter C>C used 
in one-dimensional minimization of the modified Ney/ton method, 

X ({/«,) = X^j^ -p(vA ^ S(X/^), is replaced by the matrix e^'^ - I. 


The algorithm based on equation (46) does not require inversion 

AT 

of the matrix A. Indeed, writing e in the form of a Taylor serxes, 
we obtain from (46) 
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X(t) = + tE(t) 7S(X^^) ; 
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( 47 ) 


the notation 


J i- ^ r ~r • • 

2! i! 


(48) 


is introduced. To compute the matrix E(t) we can use an algorithm 
based on application of the CayeleyHamilton theorem • The A matrix 

is a root of its characteristic equation: 
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( 49 ) 


-i- p, a"" t 

Pli) // *p, i"'' -t. .. p„ . 


The coefficients are obtained with the help of recurrent 
relationships: 


Pl= - Tj ; 

p£= — 1/2 


(50) 


- Vn -I + Tg + ... + pj V, + ), 


where Tj^ = tr (A ) . 


In (48) we limit ourselves to ra terms of the expansion. Let 
us <?all the respective matrix E^(t) . Since, according to (49), 

P(t) = 0, then 

t) = Q(t)P(t) + R(t) = R(t). (51) 

The polynomial R(t) is found as the residue of the division of the 
polynomial E (t) by the polynomial P(t). The algo.rithm makes 
it possible to find the exponential (48) to any degree of accuracy 
using exponents of the A matrix not higher than n. 

Conclusion. The routine described in this work has been 
used to process experiments since 1972. Since that time various 
improvements have been made resulting in virtually flawless opera- 
tion. Counting time under the routine is proportional td the cube of 
the number of varied parameters; it also strongly depends on the 
accuracy of the initial approximation of the parameters. The time 
of search for the statistical sum minimum for an essentially non- 
linear function f(X) comprising 5 parameters and 50 experimental 
points (of the sum of exponents type) with an EC 1040 machine is some 
15 to 20 minutes if the initial approximation was poorly given 
(the statistical sum has to be reduced by a factor of more than 10,000) , 
and two or three minutes if the statistical Siam of the initial approxi- 
mation differs from the minimum statistical sum by a factor of less 
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than 100. If the parameters of the model enter linearly, the minimum 
is found in two iterations of the Newton method, and the counting 
time decreases considerably in comparison with nonlinear models. 

In the near future it is planned to adjust a new version 
of the program, the idea of which is set forth in the last paragraph. 

The author expresses his great appreciation to V. M, Dmitriyev 
for his support and critical remarks when the program was being 
elaboratred. 
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APPENDIX. Text of Routine in FORTRAN-IV 


SUBSOUTIW 

PROCESS EXPERIMENTS BY MODIFIED NEWTON METHOD 
dimension X<*o> I <2o> »avd<io> »MARXU0> 

• ,XH1 <20> » 

*z(4A0) ,fXPiC4e«) ,yixpi(*uO} I 

• A<2a, ?o> # A1 <20'*|)> 

COMMON/TM/ M,ii|XPI<VlXPE 

comnon/tn/ n»maxx;j< 

C0MM0N/TN1/ NifXtJAVD 

commOn/helai/ ahi 
comnon/a/ a 

COHHON/GK/ 6R 
COMMON/IC/ IfOUN 
RlA0(l»10f) N 

IF(N.0T.0.AMl.N.iT.20> tOTO 1 
WRITl(S,1lO> 

STOP 

FIRST ACCESS f>RINT FUNCTION HEADING 
DO NOT COUNT FUNCTION VALUE*’ 

CALL «UN<0»2I ,NrXi FT> 

RlA0<1i109> IfS 
IF (IPl.lO.O) WXiTE<3i101» 

IFdPE.lO 1» WRITE<3»1C2> 

LOAD initial PARAMETER APPROXIMATION ARRAYS 

CALL tPOXNT(N,A|MARK,Nl fX1 .AVO.XND 
RlAD(li10A> iPt 
WRITECS#10»> |A| 
ieouN-0 

COUNT INITIAL SCALES AVD 

CALL STAT(1*I0> 

PR»EP|*8 0/l8RT<noAT<M-N41>/2.> 
call IAVO<*0,pR»ll1 »X1 ,AVDI 
MM»0 

CALL lNFORM(MM'Bo*XCOUNfN1 «XH1 »AVD> 

COUNT GRADIENTS AND SECOND DERIVATIVES MATRIX 

CALL FXSHERdpE.fl,. SO) 

DO 31 X»1 »N1 
DO 30 J*1 . N1 
AKI> J)«A(X> J) 

CONTXNUI 

ACCESS SYMMETRIC MATRIX A DIMENSION 

CAU ^NVERf (ipE.Hf.A) 

COUNT NEXT POINT BY NEWTON METHOD 
0If1"0. 

DO 10 1*1 ,M1 
F*0 . 

DO ♦ X»1iNl 
F*F*A(I . X)*eR<X> 

Dl8i*nXSl*FPiR<X) 

X1V<I)*X1 (X)-F 
Df3i»flltl/2. 

1*0 

COUNT STATISTICAL SUM AT PREDICTED POINT 

call static. ||> 

IFISS. LT.lO.lll) COTO 1» 

DO ?2 X*1.N1 
xi<x)wxi <x)doe. 

I«I41 
GOTO ?1 



HOPNEU 



*3 


*4 

c 

t 

15 
1 6 


C 


c 

t 

c 

€ 

c 

C 

1 ^ 
4 02 
103 

n’5 
;6 
1 07 

*- f: g 


IF<|»f«,ei 6010 12 

IMSS.6T.S0) eOTO 15 

MMw ? 

6610 -3 

CHECK END OF MINIMUM SEARCH 
RuSir-SS 

WRITEOMO**? 0*»iJ^- 
IF(DJ 81 .eT»AlS<R>) 

fii*so 

lMDUi!6T!c,^J!5!AfS<R».I.T.ERS*SORT<R1)) GOTO *0 

ir(fS.GI.SO) flOTO 15 

MM** 


S0“SS 

PR»EPS*S0/SORT<FHoAT<^"N+1) /2.) 

00 ^ A K»1 I N 1 

XH1 (lO«XHl <lC»*Xl CK>,*AVD<I<) 

AYDCiO«AVO<K)**«IRl'»’R/^*StAl <IC»K>> ) 


XI (K)so . 

CONTI P4UE 

TRANSFER TO NEXT COUNT STEP BY NEWTON METHOD 
GOTO « 

ONEDIMENSIONAL MINIMIZATION IN DIRECTION OF COVARIANCE- MATRIX 

DO 16 K« 1 »N 1 
GR<K)»X1 <x> 

A i 9 » / 

« w 2 V f 

AL 3 *i 


s?*s 


!3«SS , 

NED I MENS I ONAL M I N I M I ZAT I QN" SUBRQUTJ NE , 
:ALL 0DM(lPBWEPS.i2'S5»Al.<»AT.»»Ai-i 


cat 


MM*a 

goto vs 

PRINT OPTIMUM THEORETICAL FUNCTION TABLE 
dALl RTAT< 1 »SS) 

TJ °ff}WOT,T 66 uK .«1 ,*H 1 .AVS) 

PRINT FISHER MATRICES, 

COVARIANCE AND CORRELATION MATRICES, 

Accuracies normalized to statistical sum value at 

MINIMUM w . 

CALL RI 2 ULT<ie. 5 UN.M»N 1 '*SiA» 6 l»AliAV 0 »XHl) 

FORMAMIX///* IP |*0 NO CHECK PRINTOUT 
F 0 »MAT< 1 X/A/* IPl »1 - YES CHECK PRINTOUT 
F 0 RMATC 1 X»ie 0 X*'MR 

F0RMAT<19X» U# U) 

FORWAT<|ifc.T| 

fophami V.S 16 .T* Accuracy of determination of minimum', 
6 'TO scale of one -chi-square dispersion') 

FORMATC*' »Ei 6 » 7 / * predicted change of STATSUM', ■ 
*E 16 . 7 a' « REAL CHANGE OF STATSUM') 

format ( 16 ) 

formakix,' impermissible number 0 F%RAMETERS') 

GOTO 99 ■ 

RETURN 

snb 

ORIGINAL PAGE IS 
OF POOR QUALITY 


subroutine EHfXpf |H#Z,ESPI,¥i??PiJ 
OIMEHSION ZikftO) ,fXPf (%0.J #ViXPB (A00> 
READ< 1 i 160 > M 
RlAt)C*#10O (2<n*X«1rK> 

REAofi »< 01 > «|XP|{U ,I «1 
RIAOO iloO (VltXPld) fl»'. #Ki 
WRITE<3flC^> 

00 A l»1iM 

A«IXP|{I) 
ewVEXpecn 
IP<I.NE,0.> SOTO 5 
W*RITE( 3 MC 5 ) I 
STOP 


Cif 2 <I> 

URIT 6 CS» 106 ) i»A, 0 »C 
FORPATCIaJ 


sORMAT< 5 E'.P.f) 


porhat<//ix» f point number*' fe 


EXPERIMENTAL VALUE' 
COORDINATE » 

,o.MAT < 1 « , » . , 1 « . If i .7 


*' ERROR OF EXPERIMENT 
FORMATIIXi 'B« . U, » 2 ER 0 


RITURn 

tNO 


o 




/) 


SUitOUTlNi^fAVOCSfi»PR»Nj » fAVo* 

OIHfNSlON Xl f?0> .*V0<20> 

E»SII«10.E1* 

C»i^a. V 

DO 7 l«1 »N1 

CALL STATfO»SPU 

jJUpfiSSfSjfliis,) , 

xnAilt!vp5x>) •LT?ipiut a°To 5 

w5I!E<3mo«> I 
STOP 

AVBf X> bAVO(X)*C 
GOTO ? 

AVOC I>2aVOU)*SORT(PR/B> 

foRMAf^^Xf 'STATI STICAL SUM FOUND TO .BE . 

*' independent of *i5‘ 'parameter' ) 

RETURN 

END 

spcXxnt 

5|{gS’;*5 ISJlIiiriS^x. 

W»XTE<3.A) 

DO A I*1*N 
A»X<I> 

B«AVO<X> , 

IR(B.NE.O*> *0^0 7 
WRXTE<3»1d^ I 
STOP 

yRXTEC3,9» X.A.B 

REASi?“l0S» HHARMX) ,X*1 ,h) 

WR!TE<3, io4> 

Nl“*v „ 

X?<MSR|[n> 1> «0T0 IT 

Wi!TE<l.lft5> N1M 
Xl<Np«0.,.^ 

XMi<Ni)«X<X> 

AVO(Nl>»AVO<l> 

F0RMAT<!ei6 ?> 

format <iX//^ parameter! 'SXiin IT’IAL value of 
^parameter 

■ *5X. 'APPROX IMATE ERROR OF PARAMETE R ' ) '' 

P 0 RM AT <1 1 X I X 5 . E 1 S. 7 » 1 6X • E T 1 . 7 ^ ' 

FORMAt/«v//’ approx IMATE ERROR' , 

STATED equal TO ZERO’) 
FO»mAT<iX//' ALL FOLLOWING DATA PRINTOUT 
•’ACCORDING TO INTERNAL NUMERATION'/ 

•’CORRESPONDENCE OF PARAMETERS'// 

NoiMlTMOXriWTTI'^ '‘''EXTERNAL PARAMETERS') 

RETURN 

INO.. 




000 * 
000 ? 
ooc f 
000^ 
0005 
00C5 
0007 
000 * 
0009 
001 ^ 
0011 
001 7 
0011 

S§!5 
001 6 
0017 
001 8 

0019 

0020 
0021 
0022 
0028 
0024 
0029 
0024 

0027 

0028 

0029 

0030 
003’ 

0037 

0038 


4 

1 

3 


112 

113 

114 


FiSHIil 

SU 8 tOUTlNl„F|SHl»CXPErPil#S 0 > 

OTMINSXeN Xl <2g) #AyO(20>f OR<20^»tPi<20> 

COMPOW/TNI/Ni #X1 I AVD 

C6MM0N/A/ A 

COHMOK/8R/ 6R 

00 8, T«1 .N1 

Xl<f>*-1. 

CAU STAT<0»IN> 

X1?#)*1 


snv.fiiunsn <5N'*So 

IP(ABSCA(IrI)) .ST»PR/ 
AVB(1>*AV0<I>92* 

60*0 7 


> OOTO 4 


6R< T)*<SP-SNI/2 
SPt<X)«8P 
IPn.fO. 1> GOTO 1 

M*T-1 
DO 4 Kairll 

A(x I X>MSS-sMU> }-(SPB(I>«loi 
A(Krp»A(I*X> 

xijkI-o. 

ciSfivSl 

IPiXPE . 10 . 0 ) 

WPITE<3»li2> 

WilTEfS I 1 13> 

FOPMAT</f AVb ■ 8 IRf ("•PV' 

FOPMAT</» 6P»/1X, CBE14.4)> 
RETURN 
END 


GOTO 114 
<AVp<I> ,U 1 *Ni> 
, <6i(li »|»1 »N1> 
AVr*/i)(» ( 8 E 14 . 4 >> 


0001 

0007 

0008 

0004 

0005 
0004 

•0007 
000« 
0009 
001 0 
001 1 
<0012 
• 0013 

0014 

0015 
001 6 
0017 

. 0018 
001 9 
0020 
0021 
002? 
0023 

8814 
0026 
w02 7 
0028 


7 

8 


10 

11 


115 

116 
12 


((A<Xr J) iZ” 1 >N) f 


XNVERS 

subroutine INVEA8CIPE»N#A) 

DIMENSION A<20»2o> » V<20> 

N1*N“1 ^ ■ _ 

XF<TPE,EQ.0> r,OT« 3 
WPXTE(3.J1S) ■ 

DO 9 iN 

P*1./A<141> 

DO 4 I-7.N , 

VCI-<1 >«A<1 »II . 

ZfiN.EO.tX GOTO 9 
DO 8 I«1 iNl 

y**v<x)*p 

ACX fN)«y . 

DO 7 J*X#N1 

A<X .J)«A<!*1 »J9l»*V(J)*y 
CONTINUE 
A(N,N)»-P 
DO 11 I«1,N 
DO 10 J»I»N ^ 

Adf J)«-A(I*P 
A(J.I>«A(X»J> 

?SiJUur,n 

PORMATdXf'A after ACCESS '//CB|14.4n 

eIo 


O.0> OOTO 12 


VeWr'e** ^a'cCE^S^ J '<8eU?4> > 


POOf? quauty 


A <20/ 201 

jl 


ij 


/^B 


V . 

S' fi C s’ 


C'OCfi- 

tSCf 

tot? 

cos** 

&Clft 

0011 

S 81 I 
001<» 
0Q1 f 

coif 
0017 
001 i 

ml 

OOE1 

ooa? 
0021 
00 2A 

002 * 

0027 


12 


002 ' 

BBI? 110 


003 " 120 
003 ? 
0033 15 
003*. 


STAT 


Sy8«CJlttlK6 $TATCjUrS$> 

DIMF«StoN }((20)f Kl <20 )iAV 0{26) *. 


.. . 

i2<4'^O)»HXP6<«iO0ny|XPEU005 
COMMOy/TM/ M,2£IXPE» VEXf»g 
COMMON/TM/N#HAgXiX,, 

C 6 MMOK/TNI/NI rXl ,AVD 

COMMOA/HEtn^XH’ 

CDMMON/It/ICeUN 
Xi*o 

DO 5 I1*1»N 

If (HA rX<I1> .KE*0> GOTO 5 

iX^lfX., (K,>*AVD<Xi> 
continue , ^ , 

rEtfU.EO.ll enlT|C3#110) 
rcoUh*ic6uN+i 

00 <2 |»trM 

2f»2tx> • 

SI«10.Ei8 . 

fSIix^ai 

ii:aij;ssr' ’ 

IfCTU.EO.fl B«Tff3#l20> lffT»PE»Cl 


CONTINUE 
SS«S$/|. 

FoJiKtlTxJJsS?”'-* DEMONSTRATION OF rHEORETICAl CURVE*/ 
#5Xr»2 » EXPERIMENTAL VALUE'/ 

comparison with EXPERIMENT'/ 

{THEORY MINUS EXPERIMENT DIVIDED BY ERROR OF 

EXPERIMENT)'// . 

*1X. 

*53«vH->/ 

POINT I».6X»M"r6X.'r'f6Xr^2'»&X#'r»i6X**s»»6X»'t'7 

*VK# 

*53 Cl H-> > 

FORMATClX»*I' .18, X 
F0RMATtiX»53MN-)> 

RETURN 
ENO 


' »E11 kNi ' 


I ^»E11 ,A,* t t'i 
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II 


ia»S£l«»fiSS 3 WHB 


t « 


ni 


• 061 
#00? 
OOOi 
000*» 

0005 

0006 
• 007 

0009 

000*» 
00'’ n 
001^ 
001 ? 
001*5 
001 fc 
0015 
001 6 
0017 

0010 
001 <> 
OOEft 
002 ' 
002 '*' 
002’’ 
002t 
0025 

831$ 

0028 

0029 

003 0 

0031 

0032 

0033 

0034 

0039 

881$ 
0038 
0 03 9 

0040 
004*< 
0042 
0045 

0044 

0045 

0046 

004 / 


2’ 


1# 


15 


20 

21 




004fi 31 

0049 32 


0050 28 

„0051 29 

^0052 


odm 

SUmOuTlNl 

n;UU.«TrAL2) JOTO 2 
WiITi<3*30> l2^8sMt2fA|,3 
30T0 ?8 

JMS3;6T.*2> flOT# 3 

GOTO 1 

At1«-81.3 

IBSIoIitine statal{al*«s> by au finds SS 

CALI STATALfAtl »i1> , . 

IMlFf.lO.O WSlTE<3ai> Al1iS1 

Xf<S1,6t.8B> GOTO 9 

ALS-AL2 

AL2bAL1 

AL1»AL2-CAL3<.AL2>42, 

SSP«S1 

S5«82 

S2*S1 

GOTO 5 ... . 

A1.S1«<AL2^AK3? 

A2.S2*CAl,1-At3# 

A3*K3*<ALl*At2) 

6l«(Al24AL3)/2, 

B2»(AL1*A13>/2. 

AL» <Ai^ilii^|*B2*A3*B3^ ^ ?Ai-A2*A3> 

.AL.ALI (tAL2»AL*.S8fSl ^S|i93 

XP<ABs<SS"SSP> .LT.EP s*S«RT<iS> .A nO-Ii iUE.SSP) GOTO 2f 
IP(SS , Li.|l2> cOf® 15 
Zp(AL,LE.At.2l GOTO 18 
SS*Sf 

als«al 

GOTO 9 

si«st * 

8HS*8 

tPXAL.Lf ,At2» GOTO 20 

S1.S2 

AL1sA12 

GOTO 21 

SSaS2 

ALS«AL2 

AL2«Al 

S2BSS 

SSPbSS 

GOTO 0 

POAMATOK. 'INCORRECT INPUT INTO SUBROUTINE FOR 
JJjpNEDIMENSIONAL MINIMIZATION'/' 

•'82 •' #113.5/ 

•1X# 

*'*s 


J’S 

Ul 


.113 *5/ 


AL? *'#E13‘S/ 


Atl 


■' .El3-5> 


FeNMATtlXi SEARCH FOR CANONICAL CONFIGURATION'/ 

•1 X I 

•'ALl *1 *'.E13*9J 


FORMATtlX, • 
• 1 X # 

*' AL ■' #E13 
*lX. 

•' SS «' #Ei3 
STOP 
RETURN 
END, 


ARRAY OF POINTS'/ 


Si'ALi ■' f E13.5» ' AL2 ■ ' » EV3 ♦ 5 # ' AlJ 
»'#Ei3.5i' S2 ■'#I1S,5,» sss 


a I 


5#' 8i 


fits. 5/ 

IS1S.55 
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A # ■ 

V ^ V “ 
0002 
OOOf 
OOOfc 
000*5 
OOOfe 
000 ’* 
fioo« 
000 « 
001 O' 


0 C 1 * 
C 0 « i 
CO'% 

oe-* 
001 ? 
0010 

0010 


002 ^' 
0021 
0 022 


■-00* 
coc? 
0001 
0005^ 
000*^ 
000<1 
000 ? 
000? 
00 0 0 
00^ 
001 ■• 
001 ? 
oci '5 
0011* 
001 ? 
0C1? 
001 ? 
0 C 1 ? 
0 0 * ^ 
002 '* 
002 ' 

§SI? 

0 C 2 <* 

002 ? 

002 ^ 

8 BIS 

0029 

003 f'- 


IMFORM 

SUOrOuTimR IKpOrM(M»i»S 5# icoynt J.1 XHl t AVij) 
OIMFNSrON XHi(20> #AVO<20> 

IFtPM.iQ,0> t‘»lTFC3M00> 

IF(MM.60*lJ kl»lTf(3*101> 

UicajS:!! WnifliiSii 

DO 1* X*1»N1 
S*XH1 fl) 

B»AVD(I) 

WP11TE(3»166> 1»6»8 
CONTIWUE eiDCT imdutix 

FORMAT Uri^NPUT AFTER JtEP BY NEWTON METHOD') 
FeRMAT(iX#rNPUT AFTER ONEOIMENSIONAL MINIMIZATION') ' 
cORMAT{///^END. SEARCH FOR MINIMUM') 

FORMAT (1 X » E\o. *• ' “STATISTICAL SUM'. 

*»4AST BASE'} 

' * $ ^ 0 - 

FORMAKiX# 15 fOX.ilS .?f6Xi §18.?) 

RETURN 

END 


rejult 

Ai* A 1 1 1 ♦ J ) ^ ^ AV0 1 1 J ^AVO t J J J 
AI»AB*FL 0 AT<M-NRi>{*SS* 2 *> 

WRITE<3*i20> I»J*AB 

continue 

WRITE<3 ,127> 

DO A I«1|M 

WR 1 tI<( 3 * 12 «' XiJ'B 

CONTINUE 

WirTE<S*129> 

DO o IbI.N 

R»A < i R j f ^SQRY <ABS fA < I M > * A < J « J > J ) 
WRITE(3,i30> I»J'R 

RFTURN 

END 


0001 
500 ^ 
000 j 
coot 
000 ? 
SOOfc 
QOC? 
0008 
eO 09 


stataL 


CALL STATtO*fS> 

RETURN 

END 


r'AQL I 


.liAVi 




I 




An example of organizing a counting program with MODNEU subroutine /41 
call from the object module library. The FONPOL (i,Zi,N,X,FT) 
subroutine is used, which approximates the experim<?nt with a 
polynomial of the N-1 degree (Nj^2). The degree of the polynomial 
is N-1, the number of selected parameters is N. 


n ’iOQ 07ENEUT2 • USIKOV TEL. 73-25 

// HAusE ASSQN S7SRLB,X' ?9o' DISK 102 
// OPTION link 
// £Xfc ffoptran 


NAtNPtM 

CALL MOONiU 

STOP 

INS 


FUN 

SUtfRoUTXNt FUNCs,ZI,N*X,FT> 
DINknsXON AC20) 

CALL FUNP0LCZ,ZS,N*X,FTI 

RITURn 

END 


rUNPol 

SL'BKouTINE FUNPOL<X»II*N,AiFT> . 

DiMtNSlON At20> 

GOTO Z 

WRilg<Jf10«> Nl 
WRI l|(5f 101> 

GOTO 80 

2 FT»A<») 

CsAI 

00 i P«2,N 
FT»TT*C*K<il > 

3 CcCfZl 

100 format W/ 'EXPERIMENT APPROXIMATED BY POLYNOMIAL OF' 
J1'|DEGREE') ' 

101 Vo«MAT«1X>/ » X<1 >*X(2>*2 + X<5r*Z**2+X<A>*Z**^ + . . . ' > 

20 return 

ENO 

. . ... 
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